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Abstract 

Finite Volume Methods For Incompressible Flow 

by 

Darryl M. Whitlow 
Doctor of Philosophy in Applied Mathematics 

University of California, Davis 
Professor Jean- Jacques Chattot, Chair 

Two finite volume methods are derived and applied to the solution of problems of incom- 
pressible flow. In particular, external inviscid flows and boundary-layer flows are examined. 
The first method analyzed is a cell-centered finite volume scheme. It is shown to be for- 
mally first order accurate on equilateral triangles and used to calculate inviscid fiow over an 
airfoil. The second method is a vertex-centered least-squares method and is second order 
accurate. It's quality is investigated for several types of inviscid flow problems and to solve 
Prandtl's boundary-layer equations over a flat plate. Future improvements and extensions 
of the method are discussed. 
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Chapter 1 

Introduction 

Finite volume methods have been used extensively in the field of computational 
fluid dynamics. They have several advantages. Among these are their ability to be used 
on both structured and unstructured meshes. Usually, the methods can be derived in a 
straightfoward manner yet yield robust schemes with favorable properties of conservation 
of the fluxes in the flow field. 

The first part of this work concerns a finite volume method where a variable is piecewise 
constant on a triangle. The method yields a weak formulation of all necessary derivatives to 
achieve an accurate solution. The method is applied to the solution of Laplace's equation. 
The resulting difference equation for Laplace's equation is shown to have a truncation error 
that is in general no more than first order accurate on equilateral triangles. However, it will 
be demonstrated that the error of the actual solution can in general achieve an accuracy as 
high as second order on arbitrary triangles. 

In the second part of this work, a least-squares finite volume method is studied. There has 
been a resurgence of interest in least-squares type methods. These methods deal with the 
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construction of some functional whose minimization is a solution of the governing equations. 
Typically, these methods are derived in the finite element context though it has been done in 
a finite difference setting . Finite volume approaches to the construction of the functional 
are relatively young in the CFD community and started with some interest in 1996 [Q]. The 
methodology involves first integrating the governing equations over each element in the 
computational domain. Next, the functional is constructed by summing the squares of the 
discrete governing equations over each element. The functional is then minimized with 
respect to the unknowns. In general, this results in a nonlinear system which is solved 
by Newton's method. The minimization is a solution of the governing equations. Careful 
consideration must be taken to enforce the proper boundary conditions for external flow 
problems so that the resulting minimization gives a physically plausible solution. 
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Chapter 2 



A Piecewise Constant 



Approximation Method 



In this chapter, a finite volume method is derived for the solution of Laplace's 
equation in two dimensions. The method is applied to potential flow around a Joukowski 
airfoil in terms of the stream function. The solution is reached iteratively, using the Gauss- 
Seidel method, on an unstructured triangular grid. 

2.1 Derivation 

To derive the scheme, let us first suppose that the unknown ip is constant on each 
triangular element. The unknown is positioned at the circumcenter of the element. We 
seek the gradient in the normal direction across any two adjacent elements. Let the shaded 



region in figure 2.1 represent the control volume. Ci will be the path along the shaded 



region of the element to the left and C2 that of the shaded region of the element to the 



right. Integrating over this surface we have, 



Jn 



dxdy 



C 



rvi r ry2 

i^idy- j 'iIj2 dy + (b iIji dy - i/ji dy 

(^2 - V'i)(y2 - y\) 



(V^2 - ^i)(j/2 - y\) 



where the above hne integrals on C2 and Ci are zero. Similarly we can show 




Figure 2.1: 



V12 



{lp2 - 1pl){x2 - Xi) 



From and we get, 
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{ip2 - V'i)ni/i2 



where I12 is the distance between points 1 and 2. The consistency condition is Q must be 
chosen such that (V'i/'i2) " di2 = V'2 — V'l where di2 = (ii2ni and di2 is the distance between 
the two circumcenters. This imphes that, 



{'4>2 - '4>l)ll2dl2 



■02 - V'l 



n = ii2di2 = 2n 
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(2.4) 



Therefore, 
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(■02 - -01)111 



d 
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(2.5) 



Now, it is straightforward to obtain the stencil for the Laplacian. Consider a set of elements 




Figure 2.2: 

similar to those in figure 2.2 where the unknowns are again at the circumcenter of each 
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element. If we integrate the Laplacian over the center element numbered 1 we have, 



/ 



Alp dxdy 
Vip-ndl 

Cl23 

(^2-V'i)ni2 , , (V'3-V'i)ni3 , , (^4-^i)ni4 , 

— 1 • ni2£l2 H "1 • ni3t23 H -: • ni4i3i 

"12 "13 "14 

V'4 ^31 , V'3 ^23 , ^"2 ^12 , . ^31 , ^23 , ^12 x 

= —J — + — 1 — + — ^ ^^yi~ + i~ + i~) (2.6) 

"14 "13 "12 "14 "13 "12 

The local truncation error at the circumcenter of element 1 is defined as 



= A01 = (2.7) 



2.2 Numerical Verification of the Formal Accuracy 



If the elements in figure are equilateral triangles, a Taylor expansion about the 
circumcenter of the first element shows that equation ( |2.6|) is formally first order accurate 
(see appendix A). As long as all of the elements contain a circumcenter (all angles < f) 
the accuracy holds as is shown in [^. When the circumcenter is not available (i.e. right 
triangles) the centroid is used instead. 

Care should be exercised when verifying the accuracy of cell-centered type methods. This 
is due to the fact that an ordinary refinement of the grid does not preserve the orientation 
of the computational stencil and the centers from the previous grid may or may not be 
available. 



To clarify, observe the circumcenter of the center element in figure 2.3. With 



respect to equation (|2.6D, the stencil for the unknown at this point has undergone a 180° 



rotation in the refined grid of figure 2.4. Even worse, consider the case of right triangles 



Figure 2.3: First Grid 



Figure 2.4: Refined Grid 




Figure 2.5: First Grid 



Figure 2.6: Refined Grid 



where the unknowns are at the centroids. The centroids in figure 2.5 aren't centroids of any 



elements in figure 2.6 and in fact are along the edges of the triangles in the refined grid. 
Because of these possible problems the observed local truncation error of a cell-centered 
method in general will not vary monotonically with mesh refinements. Moreover, it can be 
expected to vary erratically yet, with a trend that will depend on the order of the method. 

A way to avoid this problem is to use an approach that is computationally efficient and 
can be used for cell-centered and cell- vertex methods alike. As outlined in rather than 
refining a mesh, create one computational patch. That is, create a triangulation with the 
least amount of elements that will allow the method to be fully implemented with respect 
to one unknown within the domain while enforcing a Dirichlet boundary condition. The 
truncation error is then computed on this grid and the grid is rescaled by a factor of ^ about 
the coordinates of the unknown. The truncation error is then recomputed and compared to 
that of the previous grid. This process is repeated until the formal accuracy of the method 
is apparent or machine zero is reached. This technique is used on the grid illustrated in 



figure 2^ and figure 2/I_ where the shaded triangle's centroid is given the coordinates (0,0). 



Here, the order of the truncation error is demonstrated for the following problem: 

= <f>xx + (t>yy = (2. 
where the following cases are examined: 

1. (j) = (x- 1)3 -3(x- 

2. (j) = + — Qx^y^ 
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Figure 2.7: Grid For Error Analysis On Right Triangles 

3- = sinh(7r) (sinh(7rx) sin(7ry) + sinh(7ry) sin(7r2;)) 

The circumcenter of an equilateral triangle is equivalent to it's circumcenter. Therefore, 
the variables are placed at the centroids for both equilateral and right triangles. For this 
analysis the error of the solution is defined as 

^0 — {(p 't'exactlcentroid (^'9) 

and dx is taken as the distance from the centroid of the unknown to the closest bordering 
centriod in the grid. 



In figure (2^) for a grid of equilateral triangles, the method is exact for case 1 (see 
Appendix A) and therefore is not included. The slope of the line for case 2 is exactly 1 and 
for case 3 it has a limiting slope of 1. Numerically, this confirms the formal accuracy of the 
method on a grid consisting entirely of equilateral triangles. 
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Figure 2.8: Errors corresponding to a test patch of equilateral triangles 



In figure ( p.OP for the grid consisting of right triangles, the method shows first order 
behavior for case 1 but shows second order accuracy for case 2. This is surprising since in 
general the method is not pointwise consistent on a grid consisting of right triangles. This 
means that the Laplacian from equation ( |2.6| ) will not be consistent when the variables 
are placed at the centroids. But, for these particular test cases a Taylor expansion of 



equation (2^) proves consistent with what is observed numerically. The third case is a 
true test because the Taylor expansion contains derivatives to all orders. Note that the 
truncation error is zeroth order. However, the graph of the error for that case is 
second order accurate. This illustrates an important point as discussed in [^. Namely, a 
lack of consistency does not preclude convergence. Second order accuracy of in general 
cannot be guaranteed on arbitrary triangulations. In Q second order accuracy is discussed 
for triangulations with certain regularities and symmetries. 
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-18 -16 



-14 -12 -10 
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Figure 2.9: Errors corresponding to a test patch of riglit triangles 

2.3 Inviscid Incompressible Flow Over An Airfoil 

Let us apply this method to flow over a symmetric airfoil at an angle of attack. 
Though the methodology applies to an arbitrary airfoil, here it is restricted to the class of 
airfoils that can be generated using the Karman-Trefftz conformal transformation (see ||6|). 
For a stream function ■0 = 'ip(x,y), the governing equation is 



(2.10) 



2.3.1 Numerical Boundary Conditions 



In the far field we enforce uniform flow on all elements with a node on the farfield 



boundary. Namely, 



ip = y cos a — X sin a + 



rin(r) 
2n 



(2.11) 
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where 



r = 47rasina 



(2.12) 



and a is the radius of the cyhnder used to create the airfoil via the Karman-Trefftz trans- 
formation. The angle of attack is denoted by a and is taken to be 15°. 



On the surface, the value of the streamfunction is a constant but is not known beforehand. 
It is part of the solution and must be calculated iteratively, effectively giving a Kutta 
condition. This is done by first taking the average value of the two elements immediately 
downstream of the trailing edge that are parallel to the camber line (see figure 2.1C| ) . This 



trailing edge region 




Figure 2.10: Trailing Edge 



average is then placed along the surface of the airfoil and updated after each iteration 
through the system of unknowns. For each element with an edge on the boundary, a ghost 
element is created by reflection about this edge. Now, the value at the ghost element is 
given a value such that the average of the element and the ghost element equals the value 
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on the surface. To clarify, imagine element 4 is a ghost element of element 1 in figure 2.11 
Then, 



where ^p* is the value on the surface. Hence, 



■04 = 2-0* — l/jl 



Now, the stencil we get from equation (|2.6|) can be used for all elements on the surface. 




Ghost Element 



Figure 2.11: 



2.4 Results 



The stream function was calculated for the flow over a symmetric Joukowski airfoil 
with a chord length of one and a twelve percent thickness ratio. Typical grids about the 



airfoil are shown in figures 2.12 and 2.13 were the computational domain is a [—5, 5] x [—5, 5] 
boxQ. The solutions were obtained using the Gauss-Seidel iterative method. The iterations 
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were terminated once the value of ^p* reached a plateau. As a test of convergence, the value 
of the circulation around the airfoil is compared to the exact value T = 0.888215341 for 3 
grids with varying numbers of elements on the surface of the airfoil. The results are shown 
in table |l[ The method captures the correct solution throughout the domain. A contour 



plot of the solution to the stream function is shown in figure 2.14 for a Joukowski airfoil 
with 173 surface elements. 



surface points 


r 


52 


0.945274938 


98 


0.905162059 


173 


0.880039787 



Table 2.1: Error table 



Figure 2.12: Unstructured Grid About A Joukowski Airfoil 



^AU grids produced with the aid of the "triangle" mesh generator 



Figure 2.13: Close-up of Leading and Trailing Edge 




-5 -4 -3 -2 -1 1 2 3 4 5 



Figure 2.14: Streamlines About A Joukowski Airfoil (173 surface elements) 



Figure 2.15: Close-up About A Joukowski Airfoil (173 surface elements) 
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Chapter 3 



A Least-Squares Finite Volume 
Method 



In light of the restrictions of the method derived in the previous chapter, we seek 
a higher order method that also allows greater flexibility in the type of triangulation used. 
A way to achieve this is by placing the unknowns of the governing equations at the vertices 
of the triangle and assume a linear variation over each element. 



3.1 Preliminary Derivations For A Scalar Equation 

As mentioned earlier, the governing equations are approximated over each element 



using a finite volume approximation. Figure 3.1 represents a typical triangular element for 



the method developed here. The variables are placed at the vertices which gives a linear 
approximation over each element. Consider an arbitrary variable (j) = 4'{x-, u) given at the 
three vertices. The gradient is given by: 
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J V(pdn = j> (t)ndl 
n an 
Since (f) varies linearly over the element, 



V^ri = j) (pfidl 





Figure 3.1: linear element 
Using the trapeziodal rule on each side, 



— z H ni H n2 



z z z 

where the normal vectors have been scaled to the length of the corresponding edge. Since 
for any triangle ni + n2 + n^ = 0, we have, 

It follows that the approximation for the gradient of in the triangle is, 

-1 ^ 



(3.1) 



j=i 
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Consequently, 



1 ^ 



21^ 

i=l 



-1 ^ 
i=l 

Note that since cp is linear on an element, V0 is just a constant vector. Now, consider the 
equation 

c • V(/) = (3.4) 



where c= (a, 6) = {f{x,y,(j)),g{x,y,(j))). Integrating ( |3.4| ) over an element 0^, we have, 

/ 



= c-V</)0^ (3.5) 

Where, 

j cdn (3.6) 

In general, it can be complicated to evaluate equation ( |3.6| ) exactly. However, for the 
equations approximated in this work, c varies at most linearly with respect to rc, y and (p. 
Henceforth, assuming c varies linearly on an element we get, 

3 

>: 

Combining the results of (^), (^) and (|3.7D we have. 



3 

^ y cdO = i^c,- (3.7) 



y c • Vi;^idf2 = i;^jC • nj = (3.. 

where ii^ is commonly called the cell residual. 
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3.2 The Discrete Functional 



The goal is to apply ( |3.8| ) on each triangle in the domain. For any given trian- 
gulation there are always more nodes than triangles and thus, we have more unknowns 
than available equations. To circumvent this inadequacy, one might try supplementing the 
resulting underdetermined system of equations with additional equations that conform to 
the physics of the problem and yield a unique solution |^]. In general, the process of deter- 
mining the correct additional equations is complicated. In this work a least-squares method 
is chosen instead. 

So, let the functional X be represented by, 

X=WrI^, (3.9) 

where the sum is over all of the triangles in the domain. This form of X is by no means 
obligatory but as mentioned in has the nice property that if c is a constant vector, then 
X is in fact twice the square of the L2 norm of Rt = c- 



j (c • V0)2 dn = ^ f (c • V0)^ dn = ^(c • V0)^ f 



(3.10) 



T 

2X 



In ||10[ the authors constructed a functional by calculating the residuals of all elements 
ajoining a particular node (see figure (|3.2| )). For such a functional X* we have. 



Af , N 



2 ^ ' ' 2 

i=l i=l Ti 
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where is the total number of nodes in the domain. Each element gives three contributions 

to the double summation of the right hand side. Therefore we get, 

1 ^ S 
^* = 2 ^ 2Z 4. = 2 E ^T^T = 3X (3.12) 

i=l Ti T 

So we find that either approach to building up the functional will give the same minimization 
though one can use his discretion in deciding which approach favors the intended method 
of minimization and numerical implementation. 



Figure 3.2: ilj of node i (shaded region) 

3.2.1 Minimization 

We seek a minimization of equation (|3.9| ) with respect to the variable ^. Let 

= ]^Rln^ (3.13) 

Then, 

J=^T, (3.14) 
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and 

On each element Qj, containing the node indexed i we have, 

Collecting all of these contributions at each node permits a steepest descent method of 
minimization for the discrete functional. However, in this work Newton's method was 
chosen as the method of minimization. Consequently, the Jacobian of the resulting nonlinear 
system of equations must be calculated. Let j index the nodes of an element flj,. Taking 
the partial derivative of equation (|3.16| ) with respect to (f>j gives 



d'^c , 1 dc ^ 5c ^ \ \ ^ 



2n 



Note that we have the symmetry property: 



(3.18) 



Newton's Method 



Let us call J the Jacobian whose entries are formed via equation ( |3.17| ). Let 



-z , dX dX dX .rp 

dcpi d(p2 d(t)N 



Now if 



!.i,(/)2,... ,(/.jvf (3.20) 
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using Newton's method, we solve (p from 
where q is the solution to the system: 

Jq = -f (3.22) 

The convergence criterion is either ||g|| (Euclidean norm) has to be within a desired toler- 
ance or that the calculation of T has reached a steady value to within a desired tolerance. 
Occasionally, it might be beneficial to under-relax on the correction of ( |3.21| ) in which case 
we have 



inew lold 



+ aq (3.23) 



where < a < 1. 



3.3 Two-Dimensional Systems Of Equations 

Extending the proposed least-squares method to a two-dimensional system of equa- 
tions is straightforward. Let u = (ni,M2, . . . ,Un)* G C M" and A = A(u),B = B(u) be 
n X n matrices allowing a two-dimensional n x n system to be written as 

Auo; + Buy = (3.24) 

Each element of u is piecewise linear over the elements in the domain. Integrating over an 
element 



J {Aur, + Buy)dn = Vu- J{A,B)dn\ 



(3.25) 
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We assume a linear variation of A, B with respect to u as in the scalar case. Consequently 
we get 

A = A(u) (3.26) 
B = B(u) (3.27) 

where u = ^{ui + U2 + U3). For the residual we have, 

j (Au^ + Buy) dn = (Au^ + Buy)^n^ = R^n^ (3.28) 

Now the discrete functional is 

I=IY.^^t^t (3-29) 



2 

T 



Let 



{u, v) = {{ui,Um) S u : m = 1, 2, 3, . . . , n} (3.30) 

Consider the variation of the functional on an element Qj, with respect to Ui. We have, 

^ = R* J] = R* g(Au^ + Bu^)^ 

dui ^ ^ dui ^ ^ dui 



R* T^u^ + A— ^ + — u„ + B- 

^ * OUi OUi OUi OUi 



T 



which indeed reduces to ( 3.16| ) when n = 1. As for the entries of the Jacobian we get 



/ / 5 A - dvLx dB - dvLy ^ * 



dvjdui \ \dvj ' dvj ' dvj ^ ' dvj 



T^U^ + A^:^ h 7^U„ + B- 

T 

dA - dux dB - duy 
— + A— ^ + —uy + B— ^ 

OUi dUi OUi dUi ^ ^ 

f d'^A dA dvLx dA dux 

^ \dvjdui ^ dvj dui dui dvj 

+ J!l_u ^dBd^^^dBduy\ \ 

dvjdui ^ dvj dui dui dvj J I ^ 



(3.32) 
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Chapter 4 



Inviscid Incompressible Flow 



This chapter concerns the numerical solution of Laplace's equation. Namely, 

A(/. = 

for some potential (j). This is the governing equation of inviscid incompressible flow Using 
(j)x = u and (py = V it can be written as the following system: 

= Ux+Vy 

(4.1) 

= Vx-Uy 

which are also known as the Cauchy-Riemann equations. The contribution to the functional 
on each element is 

= R^R^O^ = ^ {{ux + vyf + (vx - uyf) (4.2) 

At a node i of the element we have, 

dl 1 

-g^ = -^((^^ + ~ (^^ ~ Uy)nyJ (4.3) 

^ = ^(r^..n.,+n,,n,,.) (4.4) 

1 / ^ 

7^ — 7^ — = [nx-rinj. — ny rix ) (4.5) 
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dl 1 

- o(K + Vy)ny^ + {vx - Uy)nx^) (4.6) 



dvi 2 

a^j^ 1 



K.n^;^. - n^^ny^.) (4.7) 



— 7\ — = ^ [nvny+nx Ux ) (4. 



where j is indexes the nodes of the element. 



Equations (4.3)- ( |4.8| ) provide us with a complete set of equations for all of the 
inner nodes in the domain. The boundary nodes may not use these equations and will 
be detailed for each problem. As a numerical verification of the order of accuracy of the 
method a test is done on a unit square with Dirichlet boundary conditions. 



4.1 Test Case On A Unit Square 

Given a potential (j) = \e~^'^ sin(A;x) we have 

u = = e~^y cos{kx) (4.9) 
V = (j)y = -e-''y sm{kx) (4.10) 

where u and v satisfy the Cauchy-Riemann equations. Actually, the potential corresponds 
to the flow past an infinite sinusoidal wall. For this test, k = Qtt. If we are given Dirichlet 
boundary conditions for the variables they are implemented in a least-squares fashion by 
summing the squares of the conditions. So for this test problem, if i corresponds to a 
boundary node we have, 

^ _ {uj - e~^^' cos{kxi)f + {vj + e~^^' sin(fca;i))^ 
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giving 

di 



duidui 



Ui - e-^y^ cosikxi) (4.12) 
= 1 (4.13) 



-^ = Vi + e-^y^ sm(kxi) (4.14) 

OVi 



(4.15) 



-^ = 1 

dvidvi 

Now, the system of equations can be formulated using the following process: 

loop T = l,Tmax 

loop i £ 

if(i 3 dn) 

loop j £ Qj. 

update functional variation !F and Jacobian J using (^)- (|4.8|) 
end loop 

else update functional variation and Jacobian using ( 4.12| )- ( [4.15| ) 

end loop 
end loop 

Here, T^ax is the total number of triangles in the domain. 

Since we have Dirichlet boundary conditions, the resulting Jacobian is symmetric and 
positive definite. Therefore, Jq = —T is solved using the preconditioned conjugate gradient 
method. A Jacobi preconditioner is used. This preconditioner is defined as the matrix P 
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where 

^ , I = m 



Pi 



l,m 



(4.16) 

, otherwise 



The initial grid is a structured 13 x 13 right triangulation of the unit square (see figure 4.1). 
In this test, upon each grid refinement the solution converged in no more than 3 corrections 
with ll^l < tol = 10^^. Due to the boundary conditions and grid configuration, the method 
yielded a Jacobian in which the velocity components decoupled entirely (see Appendix B). 
Whether or not the grid was structured, the preconditioned conjugate gradient method 
consistently converged in at most, half as many iterations as there are unknowns in the 
system. In general this property cannot be guaranteed. The errors are computed using 

, = max|/(x,y)| (4.17) 



and 



(4.18) 



where r2„ is an element abutting the node i. For this problem and those following, Newton's 

i 

method is always initialized with a guess that depicts uniform flow. 



From figure 4.2 it is clear that second order accuracy is observed for this method. This 
test was repeated for larger values of k and unsurprisingly, the method was still observed 
to be second order accurate though the errors did increase for corresponding increases of 



k. For completeness a case is shown with A; = 67r on a fully unstructured grid. Table 4.1 
shows the errors in the velocity components. Figures ( [4.4D -( ^?7| ) look reasonable and can be 
improved with a finer grid. These results show we can expect the least-squares method to 
yield nice results for smooth flows if Dirichlet boundary conditions are enforced. 



Figure 4.1: Initial grid 




log2(dx) 



Figure 4.2: Errors for the sinusoidal wall 



Error 


Value 


II^M II2 

00 

ll^f II2 
loo 


0.005356007786650060 
0.036229598337966351 
0.004727775128293196 
0.025036120477625323 



Table 4.1: Error table (unstructured grid) 




Figure 4.4: u (numerical) 



Figure 4.5: u (analytical) 



Figure 4.6: v (numerical) Figure 4.7: v (analytical) 
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4.2 Inviscid Flow Over A Parabolic Profile 

In the last section the flow was reasonably smooth. Of course we are also interested 
in how the least-squares method handles flow with singular behavior. For a parabolic profile 
there is a singularity at the leading and trailing edge of the profile. The analytical solutions 
are (see pA[): 



- A 
a'' 

k 

Vexact — 9 



2k 
a 



(4.19) 
(4.20) 



where 



r2 



1 = y^{x±ay + y^ (4.21) 



and 



X ± a 



02,1 = arctan ( ) (4.22) 



For solutions computed here a = \ and k = — \ 



A- 

4.2.1 Boundary Conditions 

The computational domain is a rectangle with dimension [—2,2] x [0,2]. For v, 
Dirichlet conditions are specified everywhere and for n, everywhere except the bottom side. 
There, the functional for the governing equations is minimized with respect to u. In a 
least-squares formulation we have: 



Left, Right, Top 
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—L = {u- Uexact)i (4.24) 
OUi 

' 1 (4.25) 



duidui 



dli 

—^ = {V- Vexact)i (4.26) 

dvjdvj (4.27) 



Bottom 



v: 



u: 



dvi 
dvjdvj 



X, = ""^^""'^^ (4.28) 
dli 

{v - Vexact)i (4.29) 



dui 2 
1 



(4.30) 



^T=\ {{Ux + Vyf + [Vx - Uyf) (4.31) 

dl^ 1 



((u^ + Vy)nx^ - {vx - Uy)nyJ (4.32) 

{Ux^Ux^ + Hy^Hy^) (4.33) 



1 

— — = „ (rir Uy. — n^ rir ) (4.34) 

4.2.2 Results 

Again, the preconditioned conjugate gradient method is used to solve the system 



of equations. The initial grid is shown in figure [4.8|. From figure 4.9 observe that the error 
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does not tend to zero. The method appears to be zeroth order in l^o for v and u. This is 
due to the logarithmic behavior of u and the jump in v on the boundary y = 0. The I2 
norm of u is almost 2. In the I2 norm the order of the method for u seems to be 1 greater 
than that of v. This is surprising since v is known everywhere on the boundary. Note the 



excellent agreement in u and Uexact on the profile as illustrated in figure 4.1C for the finest 



grid. The remaining figures show contour plots corresponding to the finest grid. 




Figure 4.8: Initial Grid Used For The Computations 



'(eiTor,,) 



I (error J 



-7 -6 -5 -4 -3 

log,(dx) 



Figure 4.9: Velocity Distribution On The Lower Boundary 
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- U(numerical) 

- V 

- U(analytical) 




-2 



Figure 4.10: Velocity Distribution On The Lower Boundary 
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Figure 4.12: u (analytical) 




Figure 4.13: v (numerical) 
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4.3 Inviscid Flow Over A Cylinder 

In this section the quahty of the least-squares method is studied for flow over a 
cylinder with circulation. The flow is again governed by the Cauchy-Riemann equations. 
The analytical solutions for the velocity components are: 



liexact = 1 - 7r~9 ~ ~ (4.35) 

xF 2a^xy 

Vexact = ^ - (4-36) 



where F denotes the circulation. As usual, for all inner nodes in the domain equations (4.3)- 
El) apply. 



4.3.1 Boundary Conditions 

Farfield 

In the farfield the exact solutions are enforced: 



J. ^ {U - Uexactf + {v - Vexactfj ^^-^ 

dli 

—L = {u- Uexact)i (4.38) 
OUi 

' 1 (4.39) 



duidui 



„ {V- Vexact)i (4.40) 

OVi 

1 (4.41) 



dvidvi 



Cylinder's Surface 

For nodes on the surface, in addition to enforcing equations (|4.3| )- ( [4. 81) used in general for 
the inner nodes, contributions due to the boundary conditions are added. This is consistent 
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with the least-squares approach and provides the correct minimization. For inviscid flow 
the necessary boundary condition on the surface is u • n = 0. This is enforced at a node i 
on the surface in a least-squares formulation: 

= ^^"-'+""^'^' (4.42) 

(urix, + vny^)nx^ (4.43) 

in.f (4.44) 

rix^Uy^ (4.45) 

(urix, + vnyjuy^ (4.46) 



dui 

duidui 

dvidui 
dli 
dvi 



dvidvi 



K)' (4.47) 

(4.48) 



duidvi dvidui 

Here, the normal vector is actually of unit length in contrast to the convention established 
earlier. For a disc centered at the origin it is simply 

' (x.)2 + (y.)2 ^^-^'^ 

Enforcement of the tangency condition compromises the symmetry of the Jacobian. There- 
fore, GMRES along with the aforementioned Jacobi preconditioner (equation ( [4.16| )) is used 
within Newton's method as opposed to the preconditioned conjugate gradient method. 

4.3.2 Results 

Several cases were done for different values of the circulation. The quality with 
respect to the errors and the ability of the method to catch the correct stagnation angles 
on the cylinder were analyzed. Analytically, the stagnation angles are at 9 and n — 6 where 
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9 is given by: 



(4.50) 



For values T > AnaUoo there is one stagnation point above the cyhnder at the point (see 



(4.51) 



where 



2'naUr 



(4.52) 



For these calculations a = .5 and u^o = 1- Two grids were used in the calculations and are 



illustrated in figures 4.15-4.18. The farfield has a radius of 3 in both grids. 



From tables 4.2 and 4.3 it is clear that the errors drop approximately by a factor slightly 
smaller than ^ in each velocity component on the finer mesh. On both meshes, the errors 
increase as F was increased. Also, note that the errors in u are always larger than the 



errors in v for all values of F tested and on either grid. In Table 4.4 the data suggests 



the error of the location of the stagnation points on the surface of the cylinder increases as 



F increases on both grids and is illustrated in the streamline plots of figures 4.21-4.2E. A 



stagnation point was not found on either grid for F = A-Ka which is visible in figures 4.27] and 



4.2£. For F = Gvra the stagnation point found on the finest grid (see figures 4.3C and 4.31) 
is located at {x,y) = (0.0004,1.262) which compares favorably to the analytical values of 
(x,y) = (0, 1.3). This suggests that the difficulty of resolving the correct stagnation points 
is not simply dependent on the magnitude of F but also, whether or not the stagnation 
points lie on the surface of the cylinder. 



Figure 4.15: Grid 1 (1902 nodes 3568 elements) 




Figure 4.16: Close Up of Grid 1 (a = .5, 118 equi-spaced surface points) 



Figure 4.17: Grid 2 (5512 nodes 10548 elements) 



Figure 4.18: Close Up of Grid 2 (a 



= .5, 238 equi-spaced surface points) 
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r 


Error 


Value 





9 

M M Z 
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W^U II 


0.082418 




^1! 
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ll^f lloo 


0.042451 


27ra 


II II 2 
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ll^wlloo 


0.177046 




II II 2 


0.093906 




ll^f lloo 


0.121053 




W^u II 2 


0.175243 


2^^a^/3 


1 1 1 1 oo 


0.249885 




1 1 ^'^ 1 1 2 


0.160710 




1 1 1 1 oo 


0.189222 


Ana 


^?/, o 

M ",11^ 


0.198391 




1 1 ^« 1 1 oo 


0.276553 




Il^'y|l2 


0.185400 




oo 


0.214179 




11^^112 


0.287657 




k« oo 


0.376061 




Il^'y|l2 


0.277899 




loo 


0.310920 



r 


Error 


Value 





9 

M M Z 


0.028821 




II II QO 


0.039837 




9 


0.010564 




ll^f lloo 


0.021029 


27ra 


ll^li II2 


0.053179 




ll^wlloo 


0.082475 




II II 2 


0.044475 




ll^f lloo 


0.054971 




11^^ II2 


0.078290 


2^^a^/3 


II II 00 


0.113686 




1 1 ^'^ 1 1 2 


0.072028 




1 1 1 1 00 


0.084053 


Ana 


^11, 9 


0.087903 




ll^ulloo 


0.125113 




Il^i'll2 


0.082176 




W^v Woo 


0.094758 


6na 


II^M II2 


0.124608 




k« 00 


0.167751 




Il^'y|l2 


0.120150 




loo 


0.134961 



Table 4.2: Error table (Grid 1) Table 4.3: Error table (Grid 2) 



r 


(^analytical 


(^numerical 





0°,180° 


0°,180° 


2na 


30°, 150° 


27.4576°, 152.5424° (grid 1) 

28.7395°, 149.7479° (grid 2) 


2naV3 


60°, 120° 


54.9152°, 125.0847° (grid 1) 
57.4790°, 122.5210° (grid 2) 



Table 4.4: Error table (Grid 1 and 2) 




Figure 4.21: Streamlines (r = 0) 
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Figure 4.23: Streamlines (F = 27ra), Stagnation points at 27.4576° and 152.5424° 



Figure 4.24: Streamlines (r = 27ra\/3) 
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Figure 4.27: Streamlines (F = Ajra) 
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Figure 4.29: Streamlines (F = Ajra) 



Figure 4.30: Streamlines (F = Qira) 




Figure 4.31: Streamlines (F = Gira), Stagnation point at {x,y) = (0.0004,1.262) 
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4.4 Inviscid Flow Over An Airfoil 



In this section the Cauchy-Riemann equations are again used to approximate in- 
compressible inviscid flow over an airfoil. The treatment of the problem is much the same as 
in the previous section for flow over the cylinder. All that is needed are some modiflcations 
to the boundary conditions on the surface and the addition of a Kutta condition. We again 
will restrict the actual computations to airfoils that can be generated using the Karman- 
Trefftz conformal transformation since we have analytical solutions for the flow field. The 
methodology is nonetheless applicable to an arbitrary airfoil. However, since flow over an 
airfoil is dependent on the flow at the trailing edge, the Karman-Trefftz class of airfoils are 
good for testing the robustness of a numerical method because the velocity is always zero at 



the trailing edge. Specifically (refer to figure 4.32 ), Karman-Trefftz airfoil is chosen with a 

chord length of 1, a maximum thickness to chord ratio of 0.12, a trailing edge angle r = 10° 

and no camber. 
t=maximuni thickness 




c=chord 



Figure 4.32: airfoil nomenclature 
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4.4.1 Boundary Conditions 

Farfield 



In the farfield equations ( 4.37] )- p.41 ) are enforced where 



r 

Uexact = COS a + (4.53) 

r 

Vexact = sin a - (4.54) 



which are readily produced from the stream function of equation (2.11). As before, a is the 
angle of attack and T is the circulation given by, 

r = 47rasina (4.55) 

where a is the radius of the cylinder used to create the airfoil via the Karman-Trefftz 
transformation. For all airfoils in this section a = 0.273094. 
Airfoil's Surface 



For nodes on the surface, in addition to enforcing equations (|4.3| )- ( [4.8| ) we again enforce 
the tangency condition u • n = via equations ( 4. 42] )- (4.48). However, unlike the cylinder 



in general there is no analytic expression for the outward normal vector on the surface of 



an airfoil. Consider a node i on the surface of an airfoil as depicted in figure 4.33 . In 
order to approximate the outward normal vector at the node, first the normal vectors 
are calculated to the left and right of the node corresponding to the sides adjacent to the 
surface of triangles with the common vertex node i. Then the normal at node i is calculated 
from the arithmetic average of these two vectors. 
Kutta Condition 



With reference to figure 4.34, the Kutta condition is simply 

(ui - ua) • n = (4.56) 
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Figure 4.33: Normal vectors on the surface of an airfoil 

where ui and U2 are the average velocities in Oi and O2 respectively. The normal vector 
n is taken as the normal vector of the common edge of the two shaded elements in the 
figure. So, in addition to enforcing equations (^)- ( |4.8| ) on these two elements we need to 
minimize, 

((ui -U2) -n)^ 
^kutta = ^ (4-57) 

We only need to minimize at the nodes that aren't on the interface of and 0,2- Assume 
nodes i and j are such nodes of the respective elements ^li and Then, at node i we 
have, 

12) • n (4.58) 

(4.59) 
(4.60) 
(4.61) 
(4.62) 



d^kutta 




dUi 


3 ^ 


d'^^kutta 


_ 


duidui 
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d'^^kutta 




dujdui 


9 


d'^^kutta 




dvidui 


9 


d'^^kutta 




dvjdui 


9 
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Figure 4.34: Trailing Edge 



^^ = ^(ui-U2)-n (4.63) 

duidvi 9 

d '-^kutta IT'x^y 



dUjdVi 9 



dvidvi 9 



dvjdji 

We get similar expressions for and 



utta 



(4.65) 



kutta _ 



(4.67) 



4.4.2 Results 



Figures 4.35 and 4.36 depict the grid used in the computations. The airfoil is 



contained in a [—5,5] x [—5,5] box. Figures 4.37- 4.3£ show plots for the case a = 0. The 



results are in good agreement with the analytical solutions. Figures 



4.42 



show plots 
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for the case o; = 2° . The results show the method is able to handle flow over an airfoil with 
circulation however, it is surprising that the error in u is already becoming appreciable at 
a very small angle of attack. The error in v is mostly concentrated at the leading edge. 
Otherwise, it matches up well with the analytical solution. Adding more points on the 
surface of the airfoil and reflning the mesh did little to improve the solution. In the process 
of reaching these solutions, different formulations of the Kutta condition were enforced. 
Indeed, improving the solution signiflcantly will probably entail using a formulation of the 
Kutta condition that is more conducive to the least-squares method. 




Figure 4.35: Unstructured Grid About Airfoil 

Figure 4.36: Close-up: 70 surface points 
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Figure 4.38: u-profile about airfoil (a = 0) Figure 4.39: v-profile about airfoil (a = 0) 
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Chapter 5 



Boundary-Layer Flow 



5.1 Formulation and Derivations 

In this chapter, the least-squares method is appHed to a nonHnear set of governing 
equations in an effort to calculate boundary-layer flow over a flat plate. For high Reynolds 
regime it is known that the Prandtl boundary-layer equations coupled with potential flow 
can be used to approximate the flow as opposed to solving the Navier-Stokes equations. 
The equations are parabolic and simpler to solve. 

In dimensionless form, the Prandtl boundary-layer equations are: 



Ux + V.. 



'y 



= 



(5.1) 



UUx + VUy — Uy — UUx = 



(5.2) 



Uy — U = 



(5.3) 
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The dimensionless variables are related to the dimensional variables by: 

x' = Lx 

y' = Sy 
, UooSv 

U) = 

s 

where U' is the parallel component of velocity at the outer edge of the boundary layer, L is 
the characteristic length of the flate plate and 6 



/He 



5.1.1 Functional Minimization 



The functional for the boundary-layer equations is 

= \ ((^a: + '^2/)^ + i^^x + VUy - - UU^f + {Uy - w)^) ^l^ (5.4) 



where M is the number of elements in the domain. On each element, for each node i we 
have. 



dui 



( - 1 



— ((% - ^)nyi + {ux + Vy)nxi) ) 



2n 

V^J72-~("^^ + '"%~'^2/~^^^)(— ) (5-6) 



dujdui 
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dvi 
dujdvi 



dvjdvi 
dwjdvi 



1 ''^V /- 

-^^{'^x + Vy)nyi + -^{uUx + VUy -ojy- UUx 



+ {uUx 



vu,, 



Uy Ux UHxj + VUy. , 



UJy - UUx) 



402 



9 



-n. 



60, 



O, 



(5.9) 

(5.10) 
(5.11) 
(5.12) 



^ = [{uUx + VUy - a;, - UUx)^ - ) 0, (5.13) 



( XL TL 

{-^ + ^'f^xi + VUy Any. + — 1^ ) Q,^ (5.14) 



diOjduJi '9 402 
5.1.2 Boundary Conditions 

Bottom 



If the index i corresponds to a boundary node we have u = 0, v = 0. Consequently, 
instead of using equations (^)-(|5]l6) we have, 
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u=U, co=0 



u=l 
v=0 
00=0 



u=0, v=0 



Figure 5.1: computational domain and boundary conditions 



Vx=0 



co=0 



1i 


2 


(5.17 


dui 


= Ui 


(5.18 


duidui 


= 1 


(5.19 


dli 

dvi 


= Vi 


(5.20 


dvidvi 


= 1 


(5.21 
(5.22 



For uj, we enforce equations ( |5.13|) -( 5lT^ ). 



Left 
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Here we minimize the equation Xj 



^ '—!^ . lliis gives, 




t; — =Ui — \ 

OUi 


(5.23) 


duidui 


(5.24) 


dvi 


(5.25) 


d^i^ _^ 
dvidvi 


(5.26) 


dui 


(5.27) 


dujidbJi 


(5.28) 



Top 



On the top boundary for u and uj me minimize (" ^) = q whereby we get, 

dl- 

^=Ui-Ui (5.29) 

OUi 

' 1 (5.30) 



duidui 

For UJ we again enforce equations ( |5.27| )-( |5.28| ), and for v we use equations (|!9|)-(|1|) 



Right 



On the right side for u we simply enforce equations ( |5.5D -(5.8). For u and v we minimize 



62 




Figure 5.2: 33x65 points 



The Newton iteration is initialized with uniform flow across the plate combined 
with the no-slip condition. That is, u is set to 1 everywhere but on the plate, and v and lj 
are set to zero everywhere. Figure 5^ illustrates the grid used in these computations. On 
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the horizontal axis x varies from to 1, and on the vertical axis y varies from to 10. 

For the problem, U = 1 and Uoo = 1- Once a solution is computed it is compared to 
the analytical solution which is found by solving Blasius' equation: 

ff" + 2f"' = (5.35) 

where / = /(??). The boundary conditions are just, 

r/ = 0, / = 0, /' = 0; r/ = oo, /' = 1 (5.36) 

For the dimensionless form of the governing equations we have r/ = Also, 

u = f' (5.37) 
^ ivf-f) (5.38) 



2^ 



= ^f" (5.39) 



When ?/ = we define the shear stress r by, 



r = 0.33206^/" (5.40) 



Equation 5.35 is solved using 4th order Runge-Kutta via the Shooting Method for boundary 



value problems. 



5.2 Results 



From figure 5.3, we can conclude that the solution obtained is converged to within 
the prescribed tolerance. Although more iterations could be shown, further improvements 
in the corrective term q results in negligible improvements in the functional Z. Figure 



5.4| is also from a converged solution that is exactly the same as that corresponding to 
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figure 5^. But note that in figure 5^ the correction does not converge quadratically which 
should be the case for Newton's method. The order of the method has deteriorated because 
the residual has been compromised initially in order to reduce CPU time. Whereas the 
correction corresponding to figure was always calculated with a residual on the order of 
10^^ or less in norm, and the inner iterations on GMRES reset after 500 iterations, those 



corresponding to figure 5.4 were reset at 100 and terminated at 400 iterations regardless of 
the norm of the residual. This initially results in a residual on the order of 10"'^ but quickly 
dies down and is much less than 10~^ when the solution is converged. When executed on 
Linux machine equipped with a 550 Mhz dual Xeon Pentium, the result is a CPU time of 
353.06 seconds. This is a drastic improvement over a CPU time of 1570.08 seconds if we 



require a highly accurate solution of Newton's Method at every iteration. Figures 5.5 and 



5.6 are contour plots of u and v respectively. Figures 5.7 through 5.1C show there is very 



good agreement between the analytical and numerical results. 
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10^ 




1 2 3 4 5 6 7 8 



iteration 



Figure 5.3: Semilog plot of number of iterations versus \\q\\, T 
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Figure 5.4: Semilog plot of number of iterations versus \\q\\, T 



Figure 5.5: contour plot of u 




Figure 5.6: contour plot of v 




Figure 5.7: 




Figure 5.8: 




Figure 5.9: 




Figure 5.10: 
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5.3 Conclusions And Future Projects 

In this dissertation two numerical methods were derived. The accuracy of the 
methods was analyzed and the methods were applied to problems of incompressible flow. 

The first method studied was a cell-centered method. Though the order of the local 
truncation error for the method in general is at most first order, the order of accuracy of 
the solutions obtained with the method can be higher than first order. It was shown that 
when the method was applied to the solution of Laplace's equation, even in the test case 
that the method was not pointwise consistent a reasonable solution could be obtained. This 
can in part be traced to the fact that the truncation error never diverged in this case but 
approached a plateau. The method appears to be far more flexible than one would expect 
since it is derived by restricting the variables to be located at the circumcenter of the 
triangular elements. Relaxing this condition when necessary by placing the variables at the 
centriod of the triangular elements still tends to work reasonably well. One drawback to the 
method is that for potential flow, the velocities can only be calculated normal to the edges 
of the elements in a triangulation. Also, enforcing Dirichlet or mixed boundary conditions 
can be difficult. Nevertheless, the method is a robust alternative to solving potential flow 
problems and could easily be extended to handle advection-diffusion type equations. 

The second method studied was a least-squares finite volume method. The method was 
shown to be in general second order accurate. For some of the problems solved the method 
worked very well. However, for problems requiring a tangency boundary condition, there 
appears to be a need for further investigation into the best implementation of the boundary 
conditions. In particular, it can be troublesome for the method to resolve stagnation points 
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along such a boundary. It is also possible that solutions can be improved if as suggested in |^ 
the vertices of a triangulation are also considered variable and included in the minimization 
of the corresponding functional for a system of equations. This has the potential to allow 
for instance, characteristic grid alignment for hyperbolic systems of equations. Also, it 
seems reasonable to try constructing a least-squares method where the variables are at the 
mid-edge instead of the vertices. 

Ultimately, solutions to the Navier-Stokes equation with a least-squares finite volume 
method is desired. The methodology should be extendible to three-dimensional problems 
as well. 
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Appendix A 



Order of Accuracy For the 



Cell-Centered Method 



Proposition A. 0.1 On a grid composed of equilateral triangles 



fn Alb dxdv 

= = ^-^^ = 0{h) (A.1) 



Proof: 

Recall (refer to the figure below), 



- j 



Atp dxdy 



_ 1p4hl ^ '03^23 ^ ^^2^12 .^.^(^^l ^ hs ^ h2 

d\A di3 di2 di4, diz di2 

For a set of equilateral triangles this reduces to: 

= \/3(V'4 + V'3 + ^2-3V'i) (A.2) 



Let I equal the length of a side of a triangle. Define h = 6x = ^ and choose 
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Then Q,i = hi^y/S. With respect to the figure above, let the origin be the circumcenter of 
element numbered one and let ^' equal the solution to Laplace's equation over the domain 
of the triangulation. Then for the local truncation error we have, 



V3(*4 + -^3 + ^2- 3*1) 

k k 
^{x, y - fe) + - /i, ?/ + -) + *(x + /i, y + -) - 3*(x, y) 

y) - k^y + —%y - —%yy + ... 

( k 1 k"^ 

+ [^{x,y) - h^^ - -^y + -{h^-^xx - hk^xy + — 



Sh'^k 



xxy 



3h^ P 



(k 1 k"^ 

^{x, y) + /l*^ + -^y + -{h^^xx + hk^xy + -J^yy) 



l„o, 3h?k^ 3hk^ ^ k\ 

+ -^{h-^^xxx + -^^xxy + -^^xyy + — * 



yyy. 



-3^ix,y) 



k 



fe3 



h^k 



= —'^yy - — + h '^xx + — + —^xxv + ■■■ 



■yyy 



yy 



xxy 
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4 ■ 2 



yyy 



3 /4/i2\ 



/l3 



4 V 3 / 



^/3 V 



+ 



Division by h? completes the proof: 



by (I 



(A.4) 



From here it is easy to see why the method is exact for test case 1 in chapter 2 since for 
^ = (x - 1)3 - 3(x - l)y2 we have -^^^y = ^yyy = 0. 
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Appendix B 

On The Decoupling Of The 
Jacobian For The Cauchy-Riemann 
Equations 

For the least-squares method studied in this work, here it is shown that the re- 
sulting Jacobian used in the minimization process of the Cauchy-Riemann equations is very 
simple if the grid is composed entirely of congruent right triangles. This can result in a 
reduction in computer time if Dirichlet boundary conditions are enforced for a perspective 
problem. The resulting Jacobian will have a structure that can be exploited by a particular 
solver for large problems. It is due to the fact that the equations in the Jacobian correspond- 
ing to the minimization of any particular variable contains no entries due to contributions 
from the other variables. This is clarified by the following proposition: 

Proposition B.0.2 Let u, v be piecewise linear functions on a triangulation of a domain 

Q. Ifu,v satisfy the aforementioned discrete functional for the Cauchy-Riemann equations 
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on Q and the triangulation is composed entirely of congruent right triangles then at all inner 
nodes of the domain we have, 



— dvjdui 



where T indexes all elements abutting the node i. 



(B.ll 



Proof: 



Recall that for the Cauchy-Riemann equations we have, 



dvjdui 40 J, 



(B.2) 



where ei = (1, 0, 0), e2 = (0, 1, 0) and 63 = (0, 0, 1). Without loss of generality, assume that 
the triangles have sides of unit length and suppress the term . From equation (|B.2| ) the 



case i = j is trivial. Therefore, due to symmetry, we need only show equation (B.l) holds 
at one node where i ^ j for configurations of the figures below. 






Figure B.l: 1 config 



Figure B.2: 2"'^ config 



Figure B.3: 3'"'^ config. 
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1^* Configuration 
From equations ( p. ID and ( |B.2| ) we have, 



^ dvjdui 



((ei - rij) X nj) • es + ((-ei - n^) x nj) ■ (B.4) 
-2(n,- X n,) • eg (B.5) 
(B.6) 



2nd Configuration 



q2j 



^ dvjdui 



(rij X (ei - e2)) • eg + ((e2 - ei) x (-11^)) • eg (B.8) 
(rij X (ei - e2)) • eg - (rij x (ei - 62)) • eg (B.9) 
(B.IO) 



3'''^ Configuration 

= (rij X nj) ■ eg + ((-rij) x (-rij)) • eg (B.12) 

= (rij X rij) • eg - (rii x nj) • eg (B.13) 

= (B.14) 

Tliis completes the proof. Incidentally, this result holds if the triangulation is composed of 
congruent equilateral triangles. Referring to the figure below, 
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Figure B.4: Configuration of equilateral triangles 



Then just as in the previous configuration we have, 



q2j 

Twifr = ^ "i) ■ ^3 + (nfe X n^) ■ eg 



^ dvjdui 



= ("i X nj) • ^3 + ((-nj) X (-nj)) • 63 
= (rij X rij) • 63 - (rij x nj) • 63 



= 



(B.15) 
(B.16) 
(B.17) 
(B.18) 



This property can be used to trace errors in the coding of the numerical method. 



